toUse<-c("SacI","PvuII","NcoI","ScaI","HpaI","SpeI")
data<-read.table("TAIR_res.csv")
data2<-data[data[,2] %in% toUse,]
data3<-data2[data2[,1]=="Chr1",]
library(igraph)

minClusterSize<-100
minInsert<-1000
maxInsert<-25000


index<-1:nrow(data3)

rn<-index*2

si<-sort(data3[,4],index.return=TRUE)$ix

ssi<-sort(si,index.return=TRUE)$ix
reNumbered<-rn[ssi]
data3<-cbind(data3,reNumbered)

getEdges<-function(x,minInsert,maxInsert,data3){
 winPre<-data3[data3[,2] %in% x,c(4,6)]
 win<-sort(winPre[,1],index.return=TRUE)
 lc<-diff(win$x)
 good<-lc>=minInsert & lc<=maxInsert
 indices<-1:length(lc)
 c(sapply(indices[good],function(i) c(winPre[win$ix[i],2],winPre[win$ix[i+1],2]-
1)),recursive=TRUE)
}

sampleSizes<-function(minInsert,maxInsert){

start<-which(re.red %in% toUse)

combl2<-function(toUse,length){
 apply(t(unique(t(matrix(c(sapply(1:length(toUse),function(i){
  sapply((i+1):(i+length),function(j){
   if(j>length(toUse)){
    c(j-length(toUse),i)
   }else{
    c(i,j)
   }
  })
 }),recursive=TRUE),nrow=2)))),2,function(x) toUse[x])
}

findEnzymes<-function(x,doPlot=FALSE){
# print(x)
 if(length(unique(x))>1){
  toUse<-re.red[unique(x)]
  data3<-data[data[,2] %in% re.red & data[,1]=="Chr1",]
  index<-1:nrow(data3)
  rn<-index*2
  si<-sort(data3[,4],index.return=TRUE)$ix
  ssi<-sort(si,index.return=TRUE)$ix
  reNumbered<-rn[ssi]
  data3<-cbind(data3,reNumbered)
  l1<-sapply(toUse,getEdges,minInsert=minInsert,maxInsert=maxInsert,data3=data3)
  l1<-l1[sapply(l1,length)>0]
#  n2<-apply(combn(toUse,2),2,function(x) paste(x[1],x[2],sep="_"))
  l2<-apply(combl2(toUse,2),2,getEdges,minInsert=minInsert,maxInsert=maxInsert,data3=data3)
#  names(l2)<-n2

  g.all<-graph(c(l1,l2,recursive=TRUE),directed=FALSE)
  c.all<-clusters(g.all)
  
  linkedMarkersSum<-sum(c.all$csize[c.all$csize>2])
  largeCluster<-which(c.all$csize>=linkedMarkersSum/1000)-1
  markersToUse<-c.all$membership %in% largeCluster
  print(paste(paste(sort(toUse),collapse=" "),linkedMarkersSum,length(c.all$membership),sum(markersToUse),length(largeCluster),sep=" ## "))
  if(doPlot){
   plot((1:length(c.all$membership))[markersToUse],c.all$membership[markersToUse])
  }
  ifelse(length(largeCluster)>0,-sum(markersToUse)**3/length(c.all$membership)**2/length(largeCluster),0)
 }else{
  616
 }
}

Domains<- matrix(rep(c(1,length(re.red)),8),ncol=2,byrow=TRUE)


gres<-genoud(fn=findEnzymes, data.type.int=TRUE, BFGS=FALSE, P9=0, Domains=Domains, boundary.enforcement=2, MemoryMatrix=TRUE, max=FALSE, nvars=8, pop.size=100,starting.values=start)

c(linkedMarkers=linkedMarkersSum,markersInLargeClusters=sum(markersToUse),largeClusters=length(largeCluster))


#plot((1:length(c.all$membership))[markersToUse],c.all$membership[markersToUse])











